Persistence properties of a system of coagulating and annihilating random walkers. 
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We study a d-dimensional system of diffusing particles that on contact either annihilate with 
probability l/(q ~ I) or coagulate with probability (q — 2)/(g — f). In f-dimension, the system 
models the zero temperature Glauber dynamics of domain walls in the q-state Potts model. We 
calculate P(m,t), the probability that a randomly chosen lattice site contains a particle whose 
ancestors have undergone exactly (m — f ) coagulations. Using perturbative renormalization group 
analysis for d < 2, we show that, if the number of coagulations m is much less than the typical 
number M(t), then P(m,t) ~ m c/< V, with 6 = dQ + Q(Q - l/2)e + 0(e 2 ), C = (2Q- l)e + (2Q- 
1)(Q - l)(l/2 + AQ)e 2 + 0(e 3 ), where Q = Sri, e = 2 - d and A = -0.006 . . .. M(t) is shown 

to scale as M(t) ~ i^ 2-5 , where 5 = d(l - Q) + (Q - 1)(Q - l/2)e + C>(e 2 ). In two dimensions, 



we show that P(m, t) ~ ln(t) 



Q(3-2Q) 



ln(m) (2Q_1) i~ 2Q for m < i 2Q_1 . The f-dimensional results 



corresponding to e = 1 are compared with results from Monte Carlo simulations. 
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I. INTRODUCTION 

Persistence is understood as a property of an evolving 
system to "remember" its initial configuration for anoma- 
lously long times. A particular case of persistence that 
has received much attention is that of site persistence 
(see pj for a review). The site persistence probability is 
defined as the probability that the values of a dynami- 
cal variable at a given set of sites do not change up to 
time t. For instance, in a spin system this could be the 
probability that a spin at a given site does not flip up to 
time t, or in a reaction-diffusion system, the probability 
that no reaction takes place at that site up to time t. 
In many cases the site persistence probability decays at 
large times as a power- law 0- 

A natural generalization of site persistence is persis- 
tence of a pattern present in the initial configuration 
Q . An instance of pattern persistence would be the sur- 
vival of a test particle in a random environment. Exam- 
ples of the random environment include diffusing traps 
0, , reaction-diffusion systems such as A% <-> %A , 
A + B — > (a, L3|_or Aj+Aj — > Ai + j with mass dependent 
diffusion rates 0, 0, 113 j an d predators in predator-prey 
models 113 • The problem of calculating the survival 
probability of the test particle is hard, mainly because 
in the rest frame of the test particle, the motion of the 
other particles is correlated. 

The 1-dimensional Potts model has been a testing 
ground for various concepts of persistence. The site 
persistence problem mentioned above, has been exactly 
solved for the f-dimensional q-stdXe Potts model evolv- 
ing via zero temperature Glauber dynamics [f^ | . Besides 
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site persistence, several other persistence properties of 
the Potts model have been studied. Among these are the 
probability that a domain wall has never encountered an- 
other domain wall 0, IT3 . TTEj ] , and the probability that a 
domain present in the initial configuration survives up to 
time t [laL The former problem has been studied nu- 
merically 0, Il5| , by mean field approximations Q and 
perturbatively near q = 1 01 . However, the results ob- 
tained by these techniques do not approximate well the 
numerical results in the whole range of q. 

In dimensions greater than one, the dynamics of do- 
main walls in the Potts model is difficult to treat analyt- 
ically. Instead, we note that in f-dimension, the zero tem- 
perature Glauber dynamics of the g-state Potts model is 
equivalent to a system of diffusing particles that on con- 
tact either annihilate with probability l /(o — f) or coag- 
ulate with probability (q - 2)/{q - 1) (ll E 13. We 
study the persistence properties of this reaction-diffusion 
system in an arbitrary number of dimensions using the 
renormalization group method and calculate the expo- 
nents as an e-expansion. 

The question that we ask is, given this reaction- 
diffusion system, what is the fraction of particles that 
have never encountered another particle up to time i? 
More generally, what is the fraction of particles whose 
ancestors have undergone m coagulations up to time t? 
A convenient way to keep track of the history of coagu- 
lations is to assign a mass to each particle as follows. At 
time t = 0, let all particles be of mass one. Each time two 
particles coagulate, the new particle has a mass which is 
the sum of the masses of the two parent particles. It is 
clear that the particles of mass m will be those whose 
ancestors have undergone exactly (m — 1) coagulations. 

Let P(m, t) be the probability that a randomly chosen 
site at time t contains a particle of mass to. Let N(t) and 
p(t) denote the average particle density and average mass 
density respectively. Then, the probability distribution 



2 



P(m, t) is expected to have the scaling form 



(i) 



where l m — (m/ p(t)) 1 ^ d is the length scale associated 
with mass. The large time behavior of pit) and P{m, t) 
is characterized by two exponents 8 and 8. The mass 
density pit) ~ t~ s . For masses much smaller than the 
typical mass, P(m,t) decays as a power law in time as 
vrSl d t~ e . We will call 8 the persistence exponent. The 
exponent £ characterizes the small x behavior of the scal- 
ing function f(x) to be f(x) ~ x^ for x <C 1. Then, 
C = 2d(6 + 8 — d)/(d — 2(5), where we have used the fact 
that N(t) ~ t~ d / 2 in dimensions less than two [2(j ■ The 
two independent exponents and 8 are known for some 
limiting cases. 

When q = 2, the model reduces to the reaction diffu- 
sion model A + A — > 0. All particles are of mass one 
and it is known that P(l,t ) ~ £~ d / 2 for d < 2 and 
P(l,i) - ln(t)/t for d = 2 |13 ill. Hence, <5 = d/2, 
8 = d/2 and C = for q = 2. When q = oo, 
the model is equivalent to the reaction-diffusion system 

Ai+Aj -> [HE mm mm. smcemass isc ° n - 

served, 5 = 0. It has been shown that 8 = d+e/2 + 0(e 2 ), 
and C = e + 0(e 2 ) for e = 2 - d > 0. In 2-dimensions 
P(m,t) ~ ln(m)ln(i)/t [l0|. In 1-dimension, it is known 
via an exact calculation that P(m,t) ~ mt~ 3 ^ 2 [22| . 
When g » 1, fl has been calculated perturbatively to 
be = (q - 1)3V3/(2tt) + 0((q - l) 2 ) [14|. However for 
arbitrary values of q, the only known analytical result 
follows from a mean field approximation [2| . But the nu- 
merically obtained value for 8 is very different from the 
mean field value. In this paper, we address this issue by 
using the renormalization group formalism to systemati- 
cally calculate P(m, t) for arbitrary q. 

We now summarize our main results and give an out- 
line of the rest of the paper. In Sec. [H] we give a precise 
definition of the model and derive the stochastic partial 
differential equations obeyed by the mass distribution. 
In Sec. IIHI we express the exponent 8 in terms of the ex- 
ponent 8, though at a different value of q, reducing the 
number of unknown exponents to one. We show that 



8(q) = 8 



Thus, 



C(<z) 



2d 


8(q) + 8 


^-1, 


-d 


d-20 







(2) 



(3) 



In Sec. IIVI we use the technique developed in [Tjj to 
calculate the persistence exponent 8 as an e-expansion, 
where e = 2 — d > 0. We show that 



-dQ 
lid 



Q(Q--)e + 0(e 2 ), (4) 
where Q = If d — 2, the scaling form Eq. breaks 
down due to logarithmic corrections. We calculate these 



corrections to be 



P(m,t) 



ln(t )Q(3-2Q) ln ( m )(2Q-l) 



(5) 



given that t — > oo and m <C M(t), where M(t) is mass 
of a typical particle at time t. The analytical results for 
8 and 8 in I-dimension obtained by putting e = 1 are 
compared with the results from numerical simulations. 

In Sec. El we show that the coefficient of e" in Eq. (JIJ 
is a polynomial of degree 2n in the variable Q — (q — 
X)/q. This observation allows us to calculate the 2-loop 
corrections to the exponent £ to be 

C = (2Q-l)e+(2Q-l)(Q-l)(i + AQ)e 2 + 0(e 3 ), (6) 

where A = —0.006 . . .. The analytical results for £ in 1- 
dimension obtained by putting e = 1 are compared with 
the results from numerical simulations. 

Finally, we conclude with a summary and discussion 
in Sec. El 



II. MODEL AND FIELD THEORETIC 
FORMULATION 

In this section, we define the model and derive the 
stochastic partial differential equation obeyed by the 
mass distribution. Consider a d-dimensional lattice 
whose sites may be occupied by particles that possess 
a positive integral mass. Multiple occupancy of a lattice 
site is allowed. Given a certain configuration of particles 
on this lattice, the system evolves in time via the follow- 
ing microscopic moves, (i) With rate D, each particle 
hops to a nearest neighbor lattice site, (ii) With rate A c , 
two particles at the same site coagulates together to form 
a new particle whose mass is the sum of the masses of the 
two parent particles, (iii) With rate A a , two particles at 
the same site annihilate each other. To make connection 
with the model discussed in the introduction, we have to 
choose A c = X(q — 2)/{q— 1) and A a = X/(q— 1) where 
A is a reaction rate. The limit A — > oo corresponds to 
instantaneous reactions. In dimensions d < 2 and in the 
limit of large time, the statistical properties of a finite 
reaction rate particle system were shown to be equiva- 
lent to those of a system with infinite reaction rates ^(| . 
However, from the field theoretic point of view, it is more 
convenient to work with finite reaction rates, and hence 
A will be taken to be finite in this paper. 

Starting from the master equation for the time evolu- 
tion of the system, we now derive the effective field theory 
of the model. Let {rii} denote the configuration of parti- 
cles at site i such that n^ m is the number of particles of 
mass m at site i. Let V{. . . {nt}, {nj}, . . . ; t) be the prob- 
ability of the configuration (. . . {n^}, {rij}, . . .) at time t, 
where i and j are nearest neighbors. The master equation 
describing the time evolution of V{- ■ ■ {rii}, {rij}, . . . ;t) is 
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dP{...{n i },{n j },...) 

dt 



^2( n h™ + n O,m)'P({ni},{nj}) - y^X n i,m + l)^({"-i,m + l},{ n j,m ~ l }) 



(ij) 



^(%\m + l)V{{n i>m - 1}, {n jirn + 1}) -A c ^ ^ 



2 L m^m' 



■ /Xni,m + l)(«i,m' + ^^({^i.m + l,«i,m' + 1, n;, m+m '-l}) - ^(ni, m + 2)(nj, m + l)7 : '({n^ m + 2,n^ 2 m-l}) 



A a y^ ^ n iym ni !m ,T'({ni}) + ^ Ui, m (ni, m - l)7 : '({n i }) - ^(n^m + l)(ni >m / + l):P({ni, m + 1, n iim ' + 1}) 



EK, m + 2)(n iirra + l)P({ni, m + 2}) 



(7) 



where the time dependence of V has been dropped for 
notational simplicity and \n i rn + 1} denotes the configu- 
ration (nj,x, «i,2, ■ ■ • , ^i,m + l, ■ ■ •) at site i. The first term 
in the right hand side of Eq. Q) describes the loss and 
gain terms arising from particles diffusing to their near- 
est neighbors with rate D. The second term describes 
the loss and gain terms due to the coagulation of a pair 
of particles at a site with rate A c to form a new particle 
whose mass is the sum of the constituents. The third 
term describes the loss and gain terms due to annihila- 
tion of a pair of particles at a site with rate A a . 

The field theory corresponding to the problem can be 
derived from the master equation using Doi's formal- 
ism pfij . In short, regarding the master equation as 
a Schroedinger equation in imaginary time, the func- 
tional integral representation of the corresponding non- 
Hermitian evolution operator is constructed. This allows 
one to write down a functional integral expression for any 
correlation function of the problem, including P(m,t). 
After taking the continuum limit, one is left with the 
problem of solving an interacting field theory. The ap- 
plication of Hubbard-Stratonovich transformation to this 
field theory leaves one with a stochastic par tial differen- 
tial equation. We refer to Refs. [2^, I2H l28j for reviews 
of this procedure. Following this procedure, solving the 
master equation Eq. is equivalent to solving the fol- 
lowing Langevin equation for a stochastic field P(x, to, £): 

(d t - DV 2 )P(x,m,t) = 

- 2(A C + X a )P(x,m,t) / dm'P(x,m',t) + X C P * P 
Jo 

x,m,t), (8) 

where P * P = L dm'P(x,m' ,t)P(x,m — m',t), £ is 
white noise in space and time with unit standard devia- 
tion and i 2 = —I. 

The stochastic field P(x, to, t) is complex and is dif- 
ferent from the local mass distribution P(x,m,t), which 
denotes the number of particles of mass to in the volume 



d xdm at time t. However, the moments of P are re- 
lated to the moments of P (for instance, see [2^, l2*9jl. 

For example, P(x,m,t) = P(x,m,t), P(x,m,t) 2 — 

P(x,m,t)(Am(Ax) d )~ 1 + P(x,m,t) 2 , and so on, where 
the bar on top denotes an averaging over noise, and Aa; 
and Am, are lattice cutoffs. In this paper we only study 
the first moment of P(x,m,t), and hence disregard the 
difference between P(x, to, t) and P(x, to, t) in the rest 
of the paper. 

We will be studying the behavior of the following three 
quantities, 

P(m,x,t) for m<M(i), (9) 

p CO 

N(x,t) = / dmP(m,x,t), (10) 
Jo 

/>oo 

p(x,t) = / dm mP(m,x,t), (11) 
Jo 

where N(x, t) is the local particle density, p(x, t) is the 
local mass density and M(t) ~ p(t)/N(t) is the typical 
mass at time t. The time evolution equations obeyed by 
N(x,t) and p(x,t) are easily obtained from Eq. (JSJ. As 
for P(m, x, t), for m <C M(t), we neglect the convolution 
term in the right hand side of Eq. (JSJ) . This approxima- 
tion is justified because at large times the probability of 
collision of two light particles is negligible compared to 
the probability of collision of light and heavy particles. 
The resulting equations are: 

(d t - DV 2 )N(x, t) = -(Ac + 2A tt )iV(iT, t) 2 

+ t^2(X a + X c )^(x, t)N(x,t), (12) 

(d t - DW 2 )P(x, to, t) = -2(A C + X a )P(x, to, t)N(x, t) 
+ t^2(X c + X a )£(x, t)P(x, to, i),(13) 

(d t - DV 2 )p(x,t) = -2X a p(x,t)N(x,t) 

+ iy/2(X c + X a )£(x,t)p(x,t). (14) 
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Note that the dependence of P on mass is no longer gov- 
erned by Eq. l|13fl . Once the time dependence of P is 
calculated, its mass dependence can be restored using di- 
mensional analysis. In the rest of the paper, for the sake 
of notational simplicity, we omit the dependence of P on 
mass, unless there is a cause for confusion. 

Equations ((12(1 . 1(13(1 and 1)14(1 can be simplified as fol- 
lows. Let 



A„ 



g-2 

q-l' 
1 

~ 1' 



(15) 
(16) 



for some parameter A. Rescaling the local particle den- 
sity, local mass distribution and average density accord- 
ing to 



N(x, t),P{x, t),p(x,t) 
<7-l 



N(x,t),P(x,t),p(x,t)), (17) 



brings Eqs. (|12f| - l|14[) into the following form: 

-AiV 2 (f, t) 

+ iV2\£(x , ,t)N{x,t), (18) 
-2Q\P(x, t)N(x, t) 
+ iV2\£(x,t)P(x,t), (19) 
-2(1 - Q)Xp{x,t)N(x, t) 
+ iV2\£(x,t)p(x,t), (20) 



(fit - DV 2 )N{x , t) 
(d t -DV 2 )P(x,t) 
(d t -D\7 2 )p(x,t) 



where 



Q = 



(21) 



It can be shown that Eqs. (|18l) and 119(1 describe the two 

species reaction A + A^A, A + B^itf), in the limit 
when the concentration of A-particles is much greater 
than concentration of i?-particles. 



Secondly, there is a relation between the local mass dis- 
tribution of light particles P and the local mass density 
p. Under the substitution Q —* (1 — Q), or equivalently 
q/(q - 1) — * q, Eq. (JTHJ) transforms into Eq. (211. There- 
fore, if Fp(Q, xx, ti, x*2, t%, ■ ■ •) is a correlation function of 
P-fields, which is independent of initial conditions, then 
Fp(l — Q,xi,ti,x*2, t%, . . .) is the correlation function of 
the same configuration of p-fields. In particular, since 
P(t) ~ t~ e( Q\ we obtain p(t) ~ t- 6 ^~ Q \ Thus, we 
derive Eq. P|l. namely 



S(Q) = 9(1-Q). 



(23) 



We now examine Eqs. ((18(1 and ((19(1 for special values 
of the parameter Q. When Q = 0, the non-linear term 
in Eq. 1(191) vanishes and the concentration of monomers 
is conserved on average. Therefore, 9 — for Q = 0. 
When Q = 1/2, Eq. HHJ) is solved by P ~ N, where N 
is a solution of Eq. 1(18(1 . Then, from Eq. 1(22(1 . we obtain 
9 = d/2 for Q = 1/2. When Q = 1, it is known that 
9 = d + e/2 + 0(e 2 ), where e = 2-d 10]. If d = 1, then 
6> = 3/2, which is a consequence of an exact solution [22|] . 
rather than the e-expansion cited above. Collecting these 
results together, we have 



for Q = 0, 
for Q = i, 
§ + <3(e 2 ) forQ = l. 



(24) 



The question remains as to whether 9 for Q < 1/2 or 
equivalently q < 2 has any physical implication. For the 
site persistence problem in 1-dimension, it is known that 
the site persistence probability of the g-state Potts model 
maps to an Ising system with an initial magnetization 
given by 2/q — 1, evolving via zero-temperature Kawasaki 
dynamics |l8L Il9l 13^ . The latter system is defined at 
any value of q > 1. The correspondence between these 
two models also holds for the persistence probability of 
a single domain in the Potts model ^(|. It would be 
interesting to understand what quantity, if any, in the 
Ising model corresponds to the survival probability of 
domain walls in the Potts model. 



III. SCALING ANALYSIS OF STOCHASTIC 
EVOLUTION EQUATIONS. 

In this section, we obtain some exact results for the 
model. First, the scaled density of particles N(x, t) obeys 
the same equation as the particle density in the A + A — > 
A reaction. For this reaction, the den sity of particles 
decays for large times as t~ d / 2 in d < 2 j20,|2y. Thus, 

N ( t ) = c ^Y L ^ ( 22 ) 

where c is a constant depending on dimension only. This 
is a generalization of the exact 1-dimensional result [30l 



IV. PERTURB ATIVE COMPUTATION OF 
PERSISTENCE EXPONENT NEAR d = 2. 

In this section, we calculate the large time behavior of 
P(t) using the formalism of perturbative renormalization 
group. We closely follow the solution of the A t + Aj — > 
A i+ j model presented in Ref. . 

The solution to P(t) as a perturbative expansion in 
powers of A can be constructed from Eqs. ((18(1 and l(19|) 
using Feynman diagrams |33T ]. The Feynman rules for 
constructing terms of the expansion are summarized in 
Fig. ^ Diagrammatically, P and N are the sum of all 
Feynman diagrams with one outgoing P and N line re- 
spectively. Clearly, there are an infinite number of dia- 
grams contributing to P and N . These diagrams can be 



N 



N 




y 



-2r 



-2XQ 

FIG. 1: Propagators and vertices of the theory. 



grouped together according to the number of loops that 
they contain, thus giving rise to the loop expansion. Let 
e = 2 — d. The contribution from each diagram is a func- 
tion of the dimensionless terms XNgt and g(t) — Xt € ^ 2 and 
a term that gives the correct physical dimension [(A<) -1 
for N and {Xt)~ 2 for P}. A simple combinatorial argu- 
ment shows that the contribution from a diagram with n 
loops is proportional to g(t) n [29}. When e < 0, the main 
contribution to P and N comes from properly renormal- 
ized tree level diagrams (diagrams without loops) |34j . 
When e > 0, the loop expansion fails since for large times 
g(t) is no longer a small perturbation parameter. We 
therefore conclude that 2 is the upper critical dimension. 
For d<2we will use the formalism of perturbative renor- 
malization group to convert the loop expansion into an 
e-expansion and calculate scaling exponents as a series in 
e. 



A. Tree level diagrams 

Let iV m f and P m f be mean field densities given by the 
sum of contributions coming from tree diagrams with a 
single outgoing TV-line and P-line respectively. We de- 
note N m f and P m f by thick solid lines and thick dashed 
lines respectively. The integral equations satisfied by A m f 
and P m f are presented in diagrammatic form in Fig. [21 a) 
and^Jb). After differentiating with respect to time, they 
can be written in analytic form as 



d t P(t) 
d t N{t) 



-2QXP(t)N(t), 
-XN 2 (t), 



(25) 
(26) 



in which one can easily recognize the Smoluchowski rate 
equations of the model, obtained from Eqs. (|18|l and i|19|) 
by neglecting the noise terms in the right hand side. 
Equations (|2*5|) and i(2T)|) are easily solved yielding 



N ml (t) = 

PmfW - 

From Eq. J2HJ, we obtain 



N 



i + XN r 

Po 

(1 + XN a t) 2 Q • 



= 2Q. 



(27) 
(28) 

(29) 



(a) 



-X - 



(b) x = 



-X + 



-X + 



-2A.Q ' 



(c) 




+ 



FIG. 2: Diagrammatic form of mean field equations for (a) 
mean particle density N, (b) mean density of mass 1 particles 
P, (c) G™, and (d) G P J. 



where m f is the mean field answer for 8. 

In calculating loop corrections to any given order, we 
are faced with the problem of summing over infinitely 
many diagrams containing a given number of loops. This 
problem can be simplified by introducing mean field prop- 
agators which are sums of all tree diagrams with one in- 
coming line and one outgoing line. Expressed in terms 
of these mean field propagators, there are only finitely 
many diagrams with a fixed number of loops. 

Let (x£rf and G^f be mean field propagators. The 
integral equations satisfied by them are presented in dia- 
grammatic form in Figs.^c) and|2Id). The solutions to 
these equations are 



G£?(2|l) 
G^(2|l) 



A mf (t 2 ) 

Amf(il) 
Nmf(t 2 ) 

N mf (h) 



G (2|l), 



2Q 



G„(2|l) 



(30) 
(31) 



where 1 = 2 = (x 2 ,t 2 ) and Go is the Green's 

function of the linear diffusion equation. 



B. One loop diagrams 

Using the mean field propagators and densities, it is 
easy to classify all one loop diagrams contributing to 
P(t). These are shown in Fig. EI The computation of 
the corresponding Feynman integrals is straightforward. 
The contributions from one-loop diagrams in the limit 
-/Vq — > 00 are 



(a) 



32QXP t^ 2 



(87r) d / 2 (A At) 2£ 3e 2 (e + 2)' 



(32) 
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(a) 




(b) 




(c) 



FIG. 3: One loop corrections to the mean field result for P. 



(b) 



-64Q 2 AP £ £ / 2 

(87^)< ^ / 2 (7V A^) 2 Qe(e + 2)(e + 4) , 

-256QXP t e/2 
{8ir) d / 2 (N \t) 2 Qe 2 (e + 2) 2 (e + 4) ' 



(33) 
(34) 



where (a), (b) and (c) refer to the contributions from di- 
agrams in Figs. Efa), Etb) and Etc) respectively. Adding 
these one-loop contributions to the mean field answer 
Eq. H28|) . we obtain in the limit No — ► oo, 



P(t) = 



A 



32QXA 



37r )d/2 ei 2Q-e/5 



e + 6-2Q(e + 2) 



(e + 2) 2 (e + 4) 



+ 2- and higher loop corrections, (35) 
where A = P /(N X) 2Q . 

C. Renormalization Group Analysis of the Model. 

The large time asymptotic behavior of P(t) can be ob- 
tained by solving the Callan-Symanzik equation with ini- 
tial conditions given by Eq. (|3*5j) (see Ref. |23 for a re- 
view). The coefficients of Callan-Symanzik equation are 
determined by the law of renormalization of all the rele- 
vant couplings of the theory Eqs. I|f 8|) and (|19|l . Power 
counting analogous to that carried out in ^} , shows that 
there are only two relevant couplings of the theory in 
d < 2: the reaction rate A and the initial mass distribu- 
tion Pq . We will derive the one-loop renormalization law 
of the initial mass distribution by requiring that Eq. 1351) 
is non-singular in the limit e — ► if expressed in terms of 
renormalized relevant couplings. 

Let to be a reference time and go = Xt e J 2 be the dimen- 
sionless reaction rate. We choose to in such a way that 
go <C I. The mechanism of renormalization of the reac- 
tion rate in the theory is identical to that of the reaction 
A + A — > A. Physically, the renormalization of reaction 
rate is explained by the recurrent property of random 
walks. The probability of a reaction between particles at 
time t is proportional to the "bare" reaction rate, mul- 
tiplied by the probability that the reaction has not oc- 
curred before time t. In d < 2, the latter probability 



explicitly depends on time t. The law of renormalization 
of the reaction rate has been worked out in Ref. j^jj. If 
g R is the renormalized reaction rate, then it is related to 
go by the relation 



9R 



.9o 



1 + 9o/ 9* 



(36) 



where g* — 2ire + 0(e 2 ) is the nontrivial fixed point of 
the renormalization group flow in the space of effective 
coupling constants. The mass distribution P(to) can now 
be expressed in terms of the renormalized reaction rate 
9r to be 



Pa 



{Nog R t d /2 fQ 



1 - ^Q(2Q - 1) + 0{gl .) 



(37) 

The order g R term in Eq. (|37|l is singular at e = 0. To 
cancel this divergence, we have to introduce a renormal- 
ized initial mass distribution Pr: 

P R = Z(g R ,t ,e)P , (38) 

where Z(g R , to, e) is chosen such that 

P(t,g R ,P R ,t ) = Z(g R ,t ,e)P(t,\,P ,e), (39) 

0. Substituting Eq. ^ into 



is non-singular at e 
Eq. l|57|l. we obtain 



Z = 1 + ^Q(2Q - 1) + 0(g 2 R ) 



(40) 



The Callan-Symanzik equation is obtained by noting 
that P(t, A, Pq, e) does not depend on the reference time 
to- Therefore, 



(41) 



It follows from dimensional analysis that the most general 
form of the average mass distribution is 



(42) 



where $ is a dimensionless function. Using the scal- 
ing function Eq. (|4*2*|) in Eq. lf4*T|) . we obtain the Callan- 
Symanzik equation for P(t): 



d f3(g R ) d 
dt 2 dg R 



<K>- ) P(t,g R> P R ,t ) = 0, 

(43) 



where 



Pign) = -2t Q ^ 

Oto 



9r(.9r ~ ff*)e 
9* 



(44) 



7(3fl) = —dfo= 4^ 9R + 0(g R )m 

are the beta and gamma functions of the theory. 
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At large times, the solutions of Eq. I)43|l are governed 
by the nontrivial fixed point g* of the beta function. 
It then follows from the Callan-Symanzik Eq. Ij4ri|) that 
P(t) ~ t~ 9 where 



= Qd- 



7(5*) 



(46) 



Thus we obtain from Eqs. l(15Jl and that 



9 = dQ + Q{Q - -)e + <3(e 2 ), e > 0. (47) 

The knowledge of 9 to the first order in e combined with 
Eqs. J2J and © allows one to calculate the exponents 5 
and f with the same precision: 

5 = d(l - Q) + (Q - i)(Q - l)e + 0(e 2 ), (48) 

,2\ 



C = (2Q-l)e + 0(e z ). 



(49) 



The exponent £ is proportional to the sum of the anoma- 
lous dimensions of P and p. As a result, the mass depen- 
dence of P(to, t) can be captured neither by mean field 
theory nor by Smoluchowski approximation (see [Tof for 
a more detailed discussion of this point). 

The results of this subsection can be summarized as 
follows. The mean mass distribution P(m,t) varies as 



P(m,t) 



l/d\(2Q-l)e+0(e 2 ) 



t dQ+Q(Q-l/2)e+0(e 2 ) ■ 



(50) 



for m <C M(t), where M(t) is the typical mass at time t, 
or cquivalently the typical number of coagulations under- 
gone by all ancestors of survived particles. P(to, t) decays 
algebraically with time with an exponent independent of 
to. The coefficient multiplying this time dependent term 
does however grow algebraically with to. 



D. Two dimensions. 

The upper critical dimension of our model is 2. The 
non-trivial fixed point of the /3-function Eq. i|44|) vanishes 
at d — * 2. We therefore expect the mean field answers 
Eqs. (|27H and l|28|) to give the correct large time-small 
mass of average densities in two dimension, modulo loga- 
rithmic corrections. In this subsection we calculate these 
corrections. 

When Q = 1 it was shown that in two dimensions, 
P(m,t) ~ ln(<)ln(TO)i- 2 for t -> 00, to < M(t) [hJ. 
To calculate these corrections for arbitrary Q we need to 
solve the Callan-Symanzik equation 143|) with coefficients 
calculated at d — 2. In two dimensions, 



P{g)U=2 = 1^, 



7(flOU=2 



-2Q(2Q-l)g 
An 



0(g 2 ). 



(51) 
(52) 



Then Eq. g3| reduces to 

which has to be solved with the initial condition 



^0) = 



const 
(.9^o) 2 «' 



(54) 



provided by the mean field theory. The solution to 
Eq. (|53|) with this initial condition is 



P(i) = const x 



(ln(tAo)) Q(3 ~ 2Q) 
gf^-'hlQ 



l + 0( 



Ht/to) 



_ (55) 

When Q = 1, we recover the result of Ref. hfj. When 
Q = 0, P(t) ceases to depend on time, as it should. When 
Q = 1/2, P ~ ln(t)/t, which coincides with the decay law 
of the concentration of particles in A + A — > reaction 

The dimensional arguments that led to Eq. © cannot 
capture the mass dependent logarithmic corrections that 
are present in 2-dimensions. Hence, we need to gener- 
alize these dimensional arguments. This is provided by 
the Callan-Symanzik equation obeyed by P(m, t) when 
considered as a function of both m and t. 

The full distribution P(m, t) cannot depend on the 
choice of reference time to, which we introduced to regu- 
larize the perturbative expansion of P(t). Therefore, 



to 



dP(m,t) 
dt 



0. 



From dimensional analysis, it follows that 



P(m,t) 



N(to)*fmN(t Q ) t 



p(to) 



-F 



p(to) ' t 



,9R 



(56) 



(57) 



The form in Eq. I|57() is different from the scaling form 
used in Eq. (|42|l because Pr has to be now expressed in 
terms of to. Substituting Eq. i|57|) into Eq. (|56|l we obtain 



, , d d 



dgR 



(d p -2dN) 



P = 0, 
(58) 

where d p = 2(1 - Q) - (1 -Q)(l+2Q)g R / (Air) +0{g%) and 
dyv = 1 — gR/(An) + 0(gj i ) are scaling dimensions of fields 

and p, which can be obtained from the corresponding 
loop expansions. 

We look for solutions of Eq. I|58|l of the form 



F = Fi[ -,g R ,t 

. to 



mN(t a ) 
P(to) 



, 9r, to 



(59) 



The time dependent function Fi obeys the Callan- 
Symanzik equation H43[) ■ Using this fact and substituting 
Eq. JSHJ into Eq. (JSSJl, we obtain 



d l + Q(g R ) a , d 

^ 4 ' Wq^T) p{9r) Wr 



ram) 



^2 = 0, (60) 



8 



where T(g R ) = {2Q - l)g R /(2ird) + 0(g 2 R ). As expected, 
when d < 2, Q > | and t — > co the solution of Eq. (|rJU)) is 
given by Eq. (JSJ}. Let d = 2. Solving Eq. JSJ for Q>\ 



with the initial condition F 2 (mo) 
mean field theory, we find that 



const, provided by 



(2Q-1) 

F 2 (m) ~ (M^ 1 ) ) ( 1 + 0(1/ ln(m/m )) 



(61) 

where mo = ^jCT is a reference point in mass space. 

Combining Eqs. I|55|l and l|61(l . wc conclude that in 
d= 2 



P(m,t) 



m (^Q(3-2Q) m ( m ) 
i2Q 



(2Q-1) 



(62) 



for mo <C m <C M(t) and i — > oo. For Q = 1 we re- 
cover the answer for the average mass distribution in the 
Ai + Aj — ► Ai + j model obtained in Ref. ^(| • This result 
has also been verified numerically ^(j- When Q = 1/2, 
P(m, t) no longer depends on mass, as expected. 



E. Comparison with results from Monte Carlo 
simulations 

In this subsection, we compare the results obtained for 
the exponents 9 and S as an e-expansion with results from 
Monte Carlo simulations in 1-dimension. Let 9\ denote 
the value of 9 obtained by truncating the e-expansion at 
order e and setting e = 1. Then, 



(63) 



In Tabled we compare this analytic expression with re- 
sults from numerical simulations (see columns 2 and 3). 
There is good agreement. 

To go beyond the expression in Eq. (|63|) and to make an 
estimate of the error arising by neglecting terms of order 
e 2 and higher, we proceed as follows. Let the corrections 
from order e 2 and higher orders be denoted by i?(e, Q), 
such that 

9 = dQ + Q{Q- he + R{e,Q). (64) 

R(e, Q) must vanish at Q — and Q = 1/2 (sec Eq. (j2U)- 
Moreover, R(l, 1) = 0, since 9 = 3/2 when Q — 1 and 
e = 1 22]. Therefore, 

R(e, Q) = e 2 Q{Q- 1 -) [(1 - Q)h x {e, Q) + (e- l)h 2 (e, Q)} , 

(65) 

where h\ and h 2 are unknown functions. Setting e = 1, 
we obtain 

0=®+Q 2 + Q(Q-\)(Q-^)hi(i,Q)i (66) 



TABLE I: The numerically obtained values of 9 for different 
values of q are compared with 9\ (Eq. 1631 ~) and 62 (Eq. 1681 ). 
For q < 2, the numerical values are obtained by measuring 
the decay of mean density and then using Eq. . For q > 2, 
the numerical results are from Refs. Q and |15|. 
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Numerical 






111 
1.11 


n no ^ n m 
U.U8 ± 0.U1 


0.U6 


n no 
U.08 


1.2(3 


U.18 ± 0.U1 


A 1 /I 

0.14 


U.l < 


JL.5U 


U.6Z ± 0.U1 


0.28 


n on 
U.30 


1.77 


0.41 ±0.01 


0.41 


0.42 


2.00 


0.50 


0.50 


0.50 


3.00 


0.73 ±0.01 


0.78 


0.75 


4.00 


0.87 ±0.01 


0.94 


0.91 


5.00 


0.96 ±0.01 


1.04 


1.01 


6.00 


1.04 ±0.01 


1.11 


1.08 


8.00 


1.12 ±0.01 


1.20 


1.17 


16.00 


1.28 ±0.01 


1.35 


1.33 


25.00 


1.35 ±0.01 


1.40 


1.39 


32.00 


1.38 ±0.01 


1.42 


1.41 


50.00 


1.42 ±0.01 


1.45 


1.44 


00 


1.50 


1.50 


1.50 



The value of function h\{\, Q) at Q — can be deter- 
mined. It was shown in Ref. 14] that 



3V3 
~2T 



Q + 0(Q 2 



(67) 



Therefore, /ii(l,0) = 3y/3/n — 1. A two loop calculation 
carried out in Sec. IVl shows that the function hi(l, Q) is 
slowly varying in the interval Q £ [0,1]. Therefore, we 
replace the function h\(l,Q) by its value at Q — and 
denote the resulting expression as 9 2 . Thus, we obtain 



Q 
2 



Q 2 + Q{Q 



Z TT 



1) 



(68) 



In Table [IJ we compare 9 2 with results from Monte Carlo 
simulations in 1-dimension. The error decreases as com- 
pared to 0\. 

The error due to dropping terms of order e 2 and higher 
can be estimated. The function h\{l,Q) in Eq. I|66() is 
or order 1. The function \Q(Q — 1/2)(Q — 1)| takes on 
a maximum value of 0.05 ... in the interval Q £ [0.5, 1]. 
Hence the absolute error is of order 0.05, which is in 
agreement with the results presented in Tabled 

We do a similar analysis for 5. Let (5i be the value 
of S obtained by truncating the series Eq. H48fl and then 
putting e = 1. Then 



Si 



5Q 
2 



To obtain 5 2 , we substitute Q 
obtain 



-Q 2 . (69) 
(1 - Q) in Eq. (EH) to 



S 2 = I - ^ ± Q 2 Q(Q \){Q 1)(^ - 1). (70) 



TABLE II: The numerically obtained values of S for different 
values of q are compared with Si (Eq. 1691 ) and 82 (Eq. I7UI 1. 



q 


Numerical 


Si 


52 


2 


0.50 


0.50 


0.50 


3 


0.31 ±0.01 


0.28 


0.30 


4 


0.22 ± 0.01 


0.19 


0.22 


5 


0.18 ±0.01 


0.14 


0.17 


8 


0.11 ± 0.01 


0.08 


0.10 


16 


0.05 ±0.01 


0.04 


0.05 


00 


0.00 


0.00 


0.00 



In Table [HJ we compare the results 5\ and S 2 with re- 
sults from numerical simulations. Very good agreement 
is seen. 



(a) K=- -K + 



(b) 




+ 8 (t 2 - y 



+ 2-loop + . 



FIG. 4: (a) Schwinger-Dyson equation for P(t). (b) Pertur- 
bative expansion of the polarization operator II. 



V. THE ANALYSIS OF TWO- AND HIGHER 
LOOP CORRECTIONS. 

A. General structure of the loop expansion. 

In this subsection, we examine the contributions from 
diagrams with 2- and more loops. It will be shown that 
the coefficient of e n in the e-expansion of 8 is a polynomial 
of degree 2n in Q. It is easier to derive the result, not 
by using the formalism of renormalization group, but by 
identifying the principal set of diagrams contributing to 
the large time limit of P(t) and deriving a simple integral 
equation satisfied by the sum of these diagrams. 

The polarization operator H(t 2 ,ti) is defined as the 
sum of all one-particle irreducible diagrams with one out- 
going and one incoming P-line, with the external prop- 
agator lines stripped off. Using the polarization oper- 
ator, we can write down the Schwinger-Dyson equation 
obeyed by P(i). Let P(t) and H(t 2 ,ti) be denoted by 
a thick dashed-dotted line and by a grey circle respec- 
tively. Then, Pit) satisfies the equation shown diagram- 
matically in Fig. Ufa). In equation form, it is 

P(t) = Pud (t) + J dt 2 t\ Q J dtxU{t 2 , ti)P(ti). 

(71) 

In 2-dimensions, P{t) ~ t~ 2Q . Let rj(t) = t 2Q P(t). In 
terms of i], Eq. (|71|l reduces to 



ft rt^ 

=ilo+l dt 2 I dti 











-2Q 



. (72) 

where 770 is a constant independent of t. Differentiating 
with respect to t, we obtain 



t 2 Qu(t,h)tf Q Tj(ti). (73) 



The expansion of II(i,ti) in terms of Fcynman dia- 
grams is shown in Fig. QJb). From the Feynman dia- 
grams in Fig. it follows that the number of dotted 
lines in any given diagram is conserved as one moves 
from the right to the left. Any diagram contributing to 
the polarization operator has a single dotted line thread- 
ing through it from right to left. Also, the only vertex 
that contributes a factor Q is the PPN vertex. There- 
fore, the Q-dependent part of a given diagram has the 

form Q# oi PPN vertices n 4 (^/^+i) 2Q , where the prod- 
uct is over all vertices involving the P-line. For a k- 
loop diagram, we can have utmost 2k vertices of the 
type PPN. Also, it was shown in Ref. 01 that n-loop 
diagrams contribute at the order e™ only. Hence, after 
coupling constant renormalization, one finds that the ex- 
pression in square brackets of Eq. (|73[) is of the form 
*~ 2 J2^Li £ n P2n{Q), where P 2n (Q) is a polynomial of de- 
gree 2n in Q, and where the factor t~ 2 has been pulled 
out to give the right dimension. 

We can now solve Eq. I|73|l perturbatively, order by 
order in e. Simple dimension counting shows that 
i 2Q IT(t,ti)tj; 2Q = t^Fitx/t), where F(t) is a dimen- 
sionless function. The previous argument shows that 
F(t) — XmLi e ™P2n(Q)- Assume a power law solution 
for 7], i.e., ij — ct~ 8p where 9 p = J2 n =i a ™ e ™ an< ^ c i s a 
constant. Then Eq. I|73|) simplifies to 



9 P = - [ drF(T)T- e ». 
Jo 

Expanding F(t) and 9 p as series in e, we obtain 

— 1 J ^_ — 1 ^ — n ^' 



(74) 



711—1 



n 2 =0 



(75) 

Solving Eq. (|75|l order by order for a n , it is easy to verify 
that the coefficient of e™ is a polynomial of degree 2n in 



Q, i.e., 

oo / In \ 

n=0 \p=0 / 

where C niP 's are some unknown constants. 

Given Eq. (|76() , it is easy to rederive the one loop cor- 
rection to 9 (Eq. H47|> ~) obtained by the renormalization 
group formalism. The three unknowns in the coefficient 
of e in Eq. I|76[) ar e obtained from the exact results in 
Eq. J23 giving C lfi = 0, Ci,i = -3/2, and Ci )2 = 1. 

B. Two-loop formula for C,(Q). 

As the mean field answer for the exponent £ is 0, it 
is desirable to evaluate order e 2 correction to Eq. (|49|l . 
This requires the knowledge of order e 2 term in 9. 

From Sec. IV Al we know that the term in e 2 is a poly- 
nomial of degree 4 in Q, i.e., 

9 = dQ + Q(Q- i)e + ( £ C 2 , fc Q fc j e 2 + 0(e 3 ). (77) 
\fe=o / 

Out of the 5 unknown C^/c's, two are fixed by the condi- 
tions that 6 = for Q = and 6 = d/2 for Q = 1/2 (see 
Eq. 121). For Q = 1, it is known that 6> = tf±e/2±0(e 2 ) 
101. We assume that the order e 2 term is absent when 
Q = 1 (see Appendix for a heuristic validation of this 
assumption). This fixes the third constant and we are 
left with 

9 = dQ+Q(Q-^)e+Q(Q-^)(Q-l){AQ+B) f 2 +0(e 3 ), 

(78) 

where A and B are constants. 

The constant A is not difficult to calculate. The contri- 
bution to 9 of order e 2 Q 4 comes from the square of 1-loop 
polarization operator and from two-loop diagrams with 
four PPN- vertices. There are only two such diagrams, 
which are shown in Fig. The numerical computation 
of corresponding Feynman integrals gives 

7T 2 

A= (—- 3) - 0.296... w -0.006, (79) 

o 

where the first term on the right hand side comes from 
the order-e 2 term in the one-loop polarization operator. 

The calculation of constant B seems an almost impos- 
sible task, as there are over twenty two-loop diagrams 
contributing to it. However, the terms proportional to B 
drop out of the the two loop expression for £. Substitut- 
ing Eq. lf?%f into Eq. © one finds that 

C = (2Q-l)e+(2Q-l)(Q-l)(i + AQ)e 2 + 0(e 3 ). (80) 

In Table 11111 we compare the one-loop expression for ( 
(Eq. iPtlfll ) and two- loop expression for ( (Eq. HOJ) with 
numerical results in 1-dimcnsion. 
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FIG. 5: The two-loop diagrams contributing to the constant 
A in Eq. ffgl . 



TABLE III: Comparison of one- loop (Eq. and two- 

loop (Eq. 1801 ') results for £ with numerical simulations in 
1-dimcnsion. 



q Numerical 1-loop 2-loop 

2 0.00 0.00 0.00 

3 0.21 ±0.11 0.33 0.28 

4 0.32 ±0.08 0.50 0.44 

5 0.44 ±0.08 0.60 0.54 
8 0.59 ±0.07 0.75 0.70 

16 0.73 ±0.06 0.88 0.85 

oo 1.00 1.00 1.00 



VI. SUMMARY AND CONCLUSIONS 

In summary, we develop a systematic method to calcu- 
late the persistence exponent 9 for a system of coagulat- 
ing and annihilating random walkers, in arbitrary dimen- 
sions. In 1-dimension, this corresponds to persistence 
probabilities of domain walls in the Potts model evolving 
via zero temperature Glauber dynamics. We establish an 
exponent relation by which the number of unknown ex- 
ponents in the problem is reduced from two to one. The 
unknown persistence exponent 9 is determined perturba- 
tively using the formalism of renormalization group. 

The persistence problem studied in this paper can be 
considered as a special case of a more general problem 
of the survival probability of a test particle with diffu- 
sion constant k times the diffusion constant of the other 
particles. In this case it is known that the persistence 
exponent 9 K (q) depends on k | 3Lll4l . While simple limit- 
ing cases have been studied [Tj, |35|, Ha, H3 via numerics, 
mean-field or pcrturbative techniques, a general under- 
standing is still lacking. It would be interesting to extend 
the formalism of perturbative renormalization group to 
calculate 9 K {q). 
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APPENDIX A: MASS DISTRIBUTION IN THE 

Ai + Aj -> A i+J MODEL 

In this appendix, we present a heuristic derivation of 



the distribution of small masses in the A; + Aj 



A 



model. This model corresponds to the Q = 1 limit of the 



model discussed in the paper. For the Ai + Aj 
model, it is known H3, E2 that for m < i d/2 , 



.4 



P(m,t) 



m { t + 0(e 2 ))/d 
t d + e /2 + C){^) 

ln(m) ln(t) 

f 

1 



in d = 1 , 

in 1 < d < 2, 

in d = 2, 
in d > 2. 



(AI) 



In this appendix, we give a heuristic argument as to why 
the terms of order e 2 and higher could be absent, as a 
result of which the expansion up to order e gives the 
exact answer. 

Putting A a = in Eq. (JSJ, we obtain that the mass 
distribution P(m, x, t) evolves according to 



DV 2 ) P = \ c P*P-2\ c NP + iy / 2\ c £P, (A2) 



0_ 

di 



where N = L dmP(m) is the density of particles and 
P * P = L dm'P(m')P(m — to'). The density obeys the 
equation: 



Let 



F(s,t) 



N = X C N 2 + iy/2\ c £N. 



dmP(m, t)e 



be the Laplace transform of P(m, t). Then, 

B(s,t) = N(t) - F(s,t), 
obeys the equation 



d_ 

at 



DV 2 B = X C B 2 



(A3) 



(A4) 



(A5) 



(A6) 



The function B(s,t) obeys the same equation as the 
density N |2ja|. However, the initial conditions at t = 
are different. If the initial density of particles N(Q) = A^o, 
then B(s,0) = N (l — e~ s ). Therefore, if the average 
particle density is N(No,t), then 



F(s,t) = N(N 0> t) - N(N (1 



(A7) 



The function F(s, t) will have the scaling form 



(A8) 



where the scaling function g{x) ~ x * for x ^> 1. On 
performing the inverse Laplace transform, we obtain 



P(m,t) 



(A9) 



t d/2(l+0) ■ 

Thus, if (f> were equal to 2/d, the expression in Eq. (|A1|) 
to is exact order e. 

To calculate 4>, we look at the behavior of the particle 
density N. It is expected to have the scaling form 



N = N h(N t d/2 ) 



(A10) 



where the scaling function h(x) behaves for large x as 

x>l, (All) 



where ip is some exponent greater than zero. 

Substituting Eqs. l|A10jl and I|A11|) into Eq. fSTjl. it is 

straightforward to verify that 



(A12) 



In 1- and 2-dimensions, it is easy enough to verify that 
ip = 2 and ip = 1 respectively, consistent with the exact 
results in Eq. IjAlfl . In other dimensions, we argue as 
follows. The density of particles in the Ai + Aj — > Ai + j 
model is inversely proportional to the area swept out by 
a random walker in time t. This area varies as (y/t + c) d , 
where c is some constant. Then, the particle density 
decays as N ~ t _d / 2 (l — const/t), such that 



(A13) 



Substituting into Eq. I|A9|I . we obtain 



P(m,t) 



M-d)/d 



td/%+1 ' 

which is the one- loop answer in Eq. I|A1|1 . 



(A14) 
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